14.2 Double Machine Learning

Chernozhukov et al. (2018): ML预测 + 去偏 = -一致的因果估计

---

本节目标

• 深入理解部分线性模型与DML的理论基础
• 掌握交叉拟合（Cross-Fitting）程序
• 理解Neyman正交性为何是DML的关键
• 使用econml和doubleml库实现DML
• 通过模拟验证DML的统计性质
• 区分ML预测与因果推断的根本差异

---

ML预测vs因果推断：根本差异

危险例子：混淆变量控制

ML思维: "加入所有变量，让模型选择"

因果推断: "控制对撞变量 → 引入偏差！"

# 数据生成
X = np.random.randn(1000)
D = np.random.randn(1000)
Y = D + X + np.random.randn(1000) * 0.5  # 真实效应 = 1

# 对撞变量
Z = X + Y + np.random.randn(1000) * 0.5

# 错误：控制Z
from sklearn.linear_model import LinearRegression

model_wrong = LinearRegression()
model_wrong.fit(np.column_stack([D, X, Z]), Y)
print(f"控制对撞变量Z: beta_D = {model_wrong.coef_[0]:.3f} (有偏!)")

# 正确：只控制X
model_right = LinearRegression()
model_right.fit(np.column_stack([D, X]), Y)
print(f"只控制混淆X: beta_D = {model_right.coef_[0]:.3f} (无偏!)")

这个例子说明：即使ML模型预测精度更高（因为Z包含了Y的信息），因果效应估计却是有偏的。

---

部分线性模型

模型设定

Robinson (1988) 的部分线性模型：

其中：

• ：因果效应参数（低维，我们关心的目标）
• ：nuisance function（高维，我们需要控制但不关心的部分）
• ：倾向得分的推广（处理变量对协变量的条件期望）

朴素方法的问题

朴素ML方法：

1. 用ML估计
2. 回归对

问题：ML的正则化偏差（regularization bias）会传递到中！

当时，正则化偏差主导，不是-一致的。

---

DML的核心：Neyman正交性

Chernozhukov et al. (2018) 的关键洞察

Neyman正交条件确保nuisance参数的估计误差对目标参数的估计只有二阶影响：

即使和的收敛速度慢于也成立！

残差化（Partialling Out）

核心步骤：

Step 1：用ML估计两个nuisance functions

Step 2：构造残差

Step 3：在残差上估计因果效应

为什么"Double"？

Double有两层含义：

1. 双重残差化：同时对和做残差化（Frisch-Waugh-Lovell定理的非参数推广）
2. 双重去偏：
   • 第一重：交叉拟合消除过拟合偏差
   • 第二重：Neyman正交性消除正则化偏差

---

交叉拟合程序

为什么需要交叉拟合？

如果用同一份数据既训练ML模型又估计因果效应，会产生过拟合偏差。交叉拟合将样本分割为折（通常），每折的ML预测使用其他折的数据训练。

完整算法

DML2算法（Chernozhukov et al. 2018推荐）：

1. 将样本随机分为折：

2. 对每折：

   • 用（除第折外的数据）训练ML模型
   • 对中的观测预测：和
   • 构造残差：，

3. 汇总所有折的残差，估计：

4. 标准误：

---

Python实现

手动实现DML

import numpy as np
import pandas as pd
from sklearn.ensemble import RandomForestRegressor, GradientBoostingRegressor
from sklearn.model_selection import KFold
from sklearn.linear_model import LassoCV
import matplotlib.pyplot as plt

np.random.seed(42)

# =================================================================
# 数据生成：部分线性模型
# =================================================================

n = 2000
p = 20  # 高维协变量

# 协变量
X = np.random.randn(n, p)

# 非线性 nuisance functions
g0 = np.sin(X[:, 0]) + X[:, 1]**2 + 0.5*X[:, 2]*X[:, 3]  # E[Y|X]
m0 = 0.5*X[:, 0] + np.cos(X[:, 1]) + 0.3*X[:, 4]           # E[D|X]

# 处理变量
D = m0 + np.random.randn(n) * 0.5

# 真实因果效应
theta0 = 2.0

# 结果变量
Y = theta0 * D + g0 + np.random.randn(n) * 0.5

print(f"真实因果效应: theta_0 = {theta0}")
print(f"样本量: {n}, 协变量维度: {p}")

# =================================================================
# 朴素OLS（有遗漏变量偏差）
# =================================================================

from sklearn.linear_model import LinearRegression

naive_ols = LinearRegression()
naive_ols.fit(D.reshape(-1, 1), Y)
print(f"\n朴素OLS (不控制X): theta_hat = {naive_ols.coef_[0]:.3f} (有偏!)")

# 控制X的OLS
full_ols = LinearRegression()
full_ols.fit(np.column_stack([D, X]), Y)
print(f"控制X的OLS: theta_hat = {full_ols.coef_[0]:.3f}")

# =================================================================
# DML手动实现
# =================================================================

def dml_plr(Y, D, X, ml_model_y, ml_model_d, n_folds=5):
    """
    Double Machine Learning for Partially Linear Regression

    参数:
    ------
    Y: (n,) 结果变量
    D: (n,) 处理变量
    X: (n, p) 协变量
    ml_model_y: ML模型用于 E[Y|X]
    ml_model_d: ML模型用于 E[D|X]
    n_folds: 交叉拟合的折数

    返回:
    ------
    theta_hat: 因果效应估计
    se: 标准误
    """
    n = len(Y)
    kf = KFold(n_splits=n_folds, shuffle=True, random_state=42)

    # 存储残差
    Y_tilde = np.zeros(n)
    D_tilde = np.zeros(n)

    for train_idx, test_idx in kf.split(X):
        # 训练集
        X_train, X_test = X[train_idx], X[test_idx]
        Y_train, Y_test = Y[train_idx], Y[test_idx]
        D_train, D_test = D[train_idx], D[test_idx]

        # 拟合 E[Y|X]
        ml_y = ml_model_y.__class__(**ml_model_y.get_params())
        ml_y.fit(X_train, Y_train)
        Y_hat = ml_y.predict(X_test)

        # 拟合 E[D|X]
        ml_d = ml_model_d.__class__(**ml_model_d.get_params())
        ml_d.fit(X_train, D_train)
        D_hat = ml_d.predict(X_test)

        # 残差
        Y_tilde[test_idx] = Y_test - Y_hat
        D_tilde[test_idx] = D_test - D_hat

    # 估计 theta
    theta_hat = np.sum(D_tilde * Y_tilde) / np.sum(D_tilde**2)

    # 标准误
    residuals = Y_tilde - theta_hat * D_tilde
    sigma2 = np.mean(residuals**2 * D_tilde**2) / (np.mean(D_tilde**2))**2
    se = np.sqrt(sigma2 / n)

    return theta_hat, se

# 使用Random Forest
ml_y = RandomForestRegressor(n_estimators=200, max_depth=5, random_state=42)
ml_d = RandomForestRegressor(n_estimators=200, max_depth=5, random_state=42)

theta_dml, se_dml = dml_plr(Y, D, X, ml_y, ml_d, n_folds=5)

print(f"\nDML估计:")
print(f"theta_hat = {theta_dml:.3f} (se = {se_dml:.3f})")
print(f"95% CI: [{theta_dml - 1.96*se_dml:.3f}, {theta_dml + 1.96*se_dml:.3f}]")
print(f"真实值: {theta0}")
print(f"包含真实值? {'是' if theta_dml - 1.96*se_dml <= theta0 <= theta_dml + 1.96*se_dml else '否'}")

# =================================================================
# 使用Gradient Boosting
# =================================================================

ml_y_gb = GradientBoostingRegressor(n_estimators=200, max_depth=3, random_state=42)
ml_d_gb = GradientBoostingRegressor(n_estimators=200, max_depth=3, random_state=42)

theta_gb, se_gb = dml_plr(Y, D, X, ml_y_gb, ml_d_gb, n_folds=5)

print(f"\nDML (Gradient Boosting):")
print(f"theta_hat = {theta_gb:.3f} (se = {se_gb:.3f})")

# =================================================================
# 结果对比
# =================================================================

print("\n" + "="*60)
print("方法对比")
print("="*60)
results = pd.DataFrame({
    '方法': ['朴素OLS', '控制X的OLS', 'DML (Random Forest)', 'DML (Gradient Boosting)'],
    '估计值': [naive_ols.coef_[0], full_ols.coef_[0], theta_dml, theta_gb],
    '标准误': ['-', '-', f'{se_dml:.3f}', f'{se_gb:.3f}'],
    '真实值': [theta0]*4
})
print(results.to_string(index=False))

使用econml库

from econml.dml import LinearDML, CausalForestDML

# LinearDML
dml_linear = LinearDML(
    model_y=RandomForestRegressor(n_estimators=200, max_depth=5),
    model_t=RandomForestRegressor(n_estimators=200, max_depth=5),
    cv=5,
    random_state=42
)

dml_linear.fit(Y, D, X=X)

# 因果效应
theta_econml = dml_linear.const_marginal_effect()
ci = dml_linear.const_marginal_effect_interval(alpha=0.05)

print(f"\n[econml] LinearDML:")
print(f"theta_hat = {theta_econml[0]:.3f}")
print(f"95% CI: [{ci[0][0]:.3f}, {ci[1][0]:.3f}]")

使用doubleml库

from doubleml import DoubleMLPLR, DoubleMLData
from sklearn.ensemble import RandomForestRegressor

# 准备数据
df_dml = pd.DataFrame(X, columns=[f'X{i}' for i in range(p)])
df_dml['Y'] = Y
df_dml['D'] = D

dml_data = DoubleMLData(
    df_dml,
    y_col='Y',
    d_cols='D',
    x_cols=[f'X{i}' for i in range(p)]
)

# DoubleML
dml_plr_obj = DoubleMLPLR(
    dml_data,
    ml_l=RandomForestRegressor(n_estimators=200, max_depth=5),
    ml_m=RandomForestRegressor(n_estimators=200, max_depth=5),
    n_folds=5,
    n_rep=1
)

dml_plr_obj.fit()
print(f"\n[doubleml] DoubleMLPLR:")
print(dml_plr_obj.summary)

---

DML的统计性质

渐近正态性

在适当的正则性条件下：

其中：

双重稳健性

DML的残差化步骤类似于双重稳健估计（Doubly Robust Estimation）：

• 即使不完美，只要足够好，仍然一致
• 即使不完美，只要足够好，仍然一致
• 只需要

---

本节小结

核心要点

1. DML解决的问题：在利用ML灵活性的同时，保持因果效应估计的-一致性
2. 两个关键思想：交叉拟合（消除过拟合偏差）+ Neyman正交性（消除正则化偏差）
3. "Double"的含义：双重残差化 + 双重去偏
4. 实现简单：核心就是"残差对残差的回归"
5. Python工具：econml的LinearDML、doubleml的DoubleMLPLR

关键公式

DML估计量：

渐近分布：

---

<< 上一节：14.1 本章介绍 | 下一节：14.3 因果发现与前沿方法 >>